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Abstract. We present a case study for the semiclassical calculation of the oscillations 
in the particle and kinetic-energy densities for the two-dimensional circular billiard. 
For this system, we can give a complete classification of all closed periodic and non- 
periodic orbits. We discuss their bifurcations under variation of the starting point r and 
derive analytical expressions for their properties such as actions, stability determinants, 
momentum mismatches and Morse indices. We present semiclassical calculations of 
the spatial density oscillations using a recently developed closed-orbit theory [Roccia 
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approximations from perturbation and bifurcation theory, and test the convergence 
of the closed-orbit sum. 
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1. Introduction 

We have recently developed [U [2] a semiclassical theory for spatial density oscillations 
in quantum systems consisting of N non-interacting fermions bound in a local mean- 
field potential V(r). Our theory makes use of Gutzwiller's semiclassical approximation 
of the single-particle Green function [3] (which is essentially the Fourier transform of 
the semiclassical Van Vleck propagator [I] but adding all extra phases resulting from 
integrations done in the stationary phase approximation). While the average parts of 
the particle and kinetic-energy densities are given by the extended Thomas- Fermi theory 
(see, e.g., [5], chapter 4), their oscillating parts could be expressed P as sums over all 
closed, in general non-periodic orbits of the corresponding classical system, starting and 
ending at the same space point r. The semiclassical evaluation of the spatial density 
oscillations requires the knowledge of all closed non-periodic orbits and their actions, 
stability determinants, momentum mismatches and Morse indices (as specified in the 
next section), which in general can only be determined numerically. The resulting 
formulae are the analogues of the semiclassical trace formulae for the level density, 
initiated by Gutzwiller [3], in terms of periodic orbits. 

As an outcome of the semiclassical theory, some "local virial theorems" connecting 
kinetic and potential energy densities at any given point r and some (integro-) differential 
equations for the particle density p(r), that previously had been proved to hold exactly 
for isotropic harmonic oscillators with closed shells [6] and for linear potentials [2] in 
D spatial dimensions, could be generalized to arbitrary local potentials. Formally, they 
hold in general only in the semiclassical limit h — > 0, corresponding to N — > oo, but in 
various model potentials they have been shown [2] to be well fulfilled even for moderate 
particle numbers N. 

The present paper is a case study for the two-dimensional circular billiard, for 
which we can classify all closed orbits and determine their properties analytically, and 
that allows us to study the convergence of the orbit sum in the semiclassical formulae 
for the spatial density oscillations. In section [2] we summarize the main results of 
the semiclassical theory. In section [3] we investigate the closed orbits of the circular 
billiard. We give a complete classification of all periodic and non-periodic orbits, in- 
cluding their bifurcations that can occur at specific distances rb from the center, and 
present analytical expressions for their properties required in the semiclassical theory. 
In section H] we discuss the regularizations that become necessary at critical points where 
the semiclassical amplitudes diverge for different reasons. These are the center r = 
(symmetry breaking), the bifurcation points rb (0 < rb < R), and the boundary r = R 
(zero length of the shortest orbit). In section [5] we present results for the spatial density 
oscillations and test, in particular, their convergence and its relation to the shell effects 
in the total energy. In the appendix, we give explicit analytical results for the properties 
of some of the shortest orbits. 
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2. Semiclassical closed-orbit theory 

In this section we present the general framework and the main results of the semiclassical 
theory for spatial density oscillations developed in [H[2]. We consider a D-dimensional 
system of N non-interacting particles with mass m, obeying Fermi-Dirac statistics and 
bound by a local potential V(r). V(r) may represent the self-consistent local mean field 
of an interacting fermion system, such as it is obtained in density functional theory, or 
just a given model potential. The discrete energy eigenvalues Ej and eigenfunctions 
are given by the stationary Schrodinger equation. The quantum-mechanical 
particle density of the system at zero temperature, including a factor two for the spin 
degeneracy, is given by 



p(r) = 2^^(r)^(r), lp(r)d u r = N, (1) 



Ej<\ 



where the Fermi energy A = X(N) is determined by the normalization of the density to 
the given particle number N (which is here taken as an even integer). For the kinetic- 
energy density we discuss two different forms 

r w = -^ 2 E^w v2 ^( r )' ( 2 ) 

E 3 <\ 

^w = |^ 2 Ei v ^( r )i 2 ' ( 3 ) 

E 3 <\ 

which upon integration both lead to the exact total kinetic energy. We also investigate 
the average of these two kinetic-energy densities: 

e(r) = i[r(r)+n(r)]. (4) 

Due to the assumed time reversal symmetry and spin degeneracy, these densities are 
interrelated by the expressions 

r(p) = £(r) - \ ^ V 2 P (r) , n(r) = £(r) + \ ^ V 2 p(r) . (5) 
We now separate all densities into a smooth and an oscillating part: 

p(r) = p(r) + 5p(r), (6) 

r(r) = r(r) + *r(r), (7) 

r 1 (r) = fi(r) + <Jri(r), (8) 

e(r) =f(r) + ^(r). (9) 

The smooth parts are for differentiable potentials V(r) given by the extended Thomas- 
Fermi (ETF model) in terms of gradients of V(r). For billiard systems, where no gradient 
expansion of the potential exists, they are just given by their constant TF values. For 
the oscillating parts we have derived [1] a semiclassical expansion in terms of classical 
orbits which shall be examined analytically here for the circular billiard with radius R. 



Closed orbits and spatial density oscillations in the circular billiard 



4 



The potential of the circular billiard is defined by 

V = for < r < i? , V = oo for r> R 



(10) 



where (r, 0) are polar coordinates. We solve the Schrodinger equation for this potential 
with the Dirichlet boundary condition ipj(r = R, <p) = 0. {j} are given by the set of the 
radial quantum numbers n — 0, 1, 2, . . . and the angular-momentum quantum numbers 
I = 0, ±1, ±2, . . . The eigenenergies and normalized wave functions are given [7] by 

E nl = z 2 nl E , E = h 2 /{2mR 2 ) , 
xl>ra(r,<i>) = caJifar/fye**. (11) 

The normalization constants are c n i = [^/iiRJi + i{z n i)\~ l , where z n \ is the n-th zero of 
the cylindrical Bessel function Ji(z). ip n i can be summed to yield the exact quantum- 
mechanical densities (CQ) - ([3]). 

The smooth parts of the densities, which in billiard systems are independent of r, 
are given by their TF values^! in terms of a smooth Fermi energy A: 



P = Ptf 



mX 



7T 



h 2 2ttR 2 




mX 2 



ttf 



2ttH 2 AttR 2 




(12) 



Note that f is the common smooth part of all three kinetic energy densities (J2J) - (HJ)- 
From the Weyl expansion [HI [9] of the integrated level density given in fl86l) , we can find 
the asymptotic expansion of the Fermi energy X(N) in powers of N" 1 ^ 2 : 



X(N) = E 



2N + 2(2N) 1/2 + - + 0(N^ 2 ] 
3 



(13) 



For the oscillating parts of the densities, we now reproduce the main formulae 
from [1] for the special case of D = 2 space dimensions and spherical symmetry when 
all densities only depend on the radial variable r = |r|. To leading order in h, the 
semiclassical expressions for the oscillating parts of the densities are given by 

1 „ ~ s 7T 37T 



Sp(r) 
5r(r) 



r cos 



h 



7 
2 

— y^AJX,r) cos 
2m ^ <y ' 



-5 7 (A,r) 



^7" 



4 



^5 7 (A,r) - /i 7 | 



3tt 
T 



2 

Sn(r) ~ ^^ 7 (A,r)g 7 (A,r) cos 



2m 



1 ~ it 3tt 

-S 7 (X,r) 



(14) 
(15) 
(16) 



The sum is over all closed orbits 7 starting and ending in the point r. For periodic 
orbits (POs) the action integral S" 7 becomes independent of r. Therefore, the POs do 
not yield any oscillating phases in the above expressions and their contributions are 
smooth functions of r varying only through the amplitude factors *4 7 and Q 7 in the 
expressions above. The leading contributions to the density oscillations therefore come 

I We note that the density ptf hi (TT^)) cannot be normalized to the correct particle number N. The 
reason is that the constant TF density is not able to reproduce the sharp decrease of the quantum- 
mechanical density p{r) near the boundary r = R, where it is forced to become zero. 
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from the non-periodic orbits (NPOs). In one-dimensional systems it has, in fact, been 
shown [H[2] that the contributions of the POs are completely absorbed by the smooth TF 
parts of the densities. In higher- dimensional systems, like the present circular billiard, 
some contributions of POs must be included in connection with symmetry breaking at 
r = and with bifurcations at finite distances r > 0, as will be discussed in section HI 

The action function S 1 (X,r) = S 1 (X,r,r' = r) is gained from the general open 
action integral for an orbit starting at r and ending at r' at fixed energy E = X 



S y (X,r,r') = J p(A,q)-dq, (17) 
where p(A, r) is the classical momentum in the point r, for V(r) = given by 



p(A,r) = ^-V2mX, (18) 



r 



whose modulus is a constant of motion and denoted here by p\ = |p(A, r)|. /i 7 is the 
Morse index that counts the number of conjugate points along the orbit [3] and will be 
discussed in section |3~T1 The quantity Q 7 (A,r) appearing in (|T6|) for 5ri(r) is defined as 

Q 7 (A, r) = [P(^)- P (A,r')] 7 ,^ r = cog[ ^ p% (ig) 

P\ 

where p and p' are the short notations for the initial and final momentum, respectively, 
of the orbit 7 at the point r. These are obtained also from the action integral fTT7|) by 

V rl S 7 (A,r,r / ) = -p , V r ^ 7 (A,r,r') = p' . (20) 

r=r' r=r' 

Since Q 7 (A,r) in (fT9|) depends on the angle 9 between p and p', it may be called the 
"momentum mismatch function", being Q 7 = +1 for POs with p = p' and Q 7 = — 1 
for self-retracing NPOs with p = — p'. 

The common semiclassical amplitude A^(X,r) in all densities is given by 



(21) 



np\ V2^T 7 (A,r) ' \dr[ 

Here T 7 (A, r) = d5' 7 (£', r)/dE\ E= ^ is the running time of the orbit, and Pj_ 7 its stability 
determinant calculated from the components p± and r[_ transverse to the orbit 7 of the 
initial momentum and final coordinate, respectively. We shall in the following express 
the amplitudes defined in ( |2T!) through the Jacobians Jly, which are defined as the inverse 
stability determinants, and omit the Fermi energy A from all arguments: 

, / x 2m 1 , , f 5r\ 

Ay(r) = , J^ r ) 



From (JT4ll and (fT5l) using p\/2m = A, we find immediately the semiclassical relation 

5r(r) ~ X6p(r), (23) 

which is the local virial theorem discussed extensively in j2] (valid here for V(r) = 0). 
Note, however, that this theorem fails at the boundary r = R (cf. section fl~3l) . 
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It has been observed in [2] (see also figures [T^] and [12]) that for systems with spherical 
symmetry for D > 1 (except harmonic oscillators), there exist two types of oscillations: 

(i) regular, short-ranged oscillations due to the orbits librating in the radial direction, 
which we will briefly call the radial orbits, and 

(ii) irregular, long-ranged oscillations due to the non-radial orbits. (Note that in 
isotropic harmonic oscillators, all non-radial orbits are periodic, which explains the 
fact [21 [6] that these systems do not exhibit any irregular oscillations.) Like in [2], we 
therefore decompose the oscillating parts of the densities ( 1141) - ( fl~6j) in the following 
way (in obvious notation) 

5p(r) = S T p(r) + 5 iTT p(r) , 
5r(r) = <5 r r(r) + 5 irr r(r) , 

5r 1 (r) = 6 I T 1 (r) + 6 in: T 1 (r). (24) 

The contribution of the radial NPOs to the particle density (114j) near r = can 
be summed over all repetitions k and has already been given in [H [2]. For the present 
system it becomes: 

8 T p(r) = (-lf- 1 -^- J (2rp x /h) . (r « R) (25) 
hT rl (X) 

Here Jq{z) is the cylindrical Bessel function of order zero, M is the number of filled "main 
shells" to be discussed in section [5], and T T \ = 4mR/p\ is the period of the diametrical 
PO. This result will be generalized to larger radii in section 14.11 The regularity of 
the short-ranged oscillations is due to the fact that all radial NPOs including their 
repetitions k (see their discussion in section [331) contribute with the same period in r at 
all distances except near the boundary r = R. This is different for the non-radial NPOs 
(which will be discussed in section 13. 6ft : their periods in r depend on their individual 
forms and become larger in the limit r — > (cf. also section |4.1. 2p . which explains the 
irregular nature of the long-ranged oscillations. 

The contribution of the radial NPOs to the kinetic energy densities are then easily 
obtained from the general equations ( fl~4l) - ( fT6l) . For 5 r r(r) we may just use the local 
virial theorem ( 1231) to get an expression valid for small r. Since Q 1 {r) = — 1 for all 
radial NPOs, we get the relation [1] 

5 r ri(r) = -5 T r(r) , (r < R) (26) 

again valid for small r. 

In the following section, we shall classify all closed orbits in the circular billiard 
and derive their analytical properties introduced above. 
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3. Closed orbits in the circular billiard 

3.1. Classification and symmetries of orbits 

We may classify the NPOs in the circular billiard in the same way as it has been done 
for the POs by Balian and Bloch [10]. They specified them by pairs of integers (v,w), 
where v is the number of reflections at the boundary and w is the winding number 
around the center (r = 0). If v and w have a common divisor k > 1, the number k is 
the number of repetitions of the primitive orbit (which by definition has k — 1). The 
librating diametrical POs have by definition the winding number v = w. For the non- 
radial NPOs, there generally exist no repetitions (i.e., k = 1), except at the points r = 
and r = R where they become periodic. Repetitions of non-radial NPOs do occur at 
specific isolated points < r < R at which they become identical with POs or fractions 
thereof. In the calculation of the semiclassical densities we found, however, that these 
repetitions are practically negligible because they lead to rather long orbits. For the 
radial NPOs, a generalized repetition number k > will be defined in section 13.61 

All orbits starting and ending in a fixed point r > are isolated and occur in discrete 
degenerate pairs corresponding to the time-reversal symmetry. This holds also for the 
librating (i.e., radial) POs|§| All orbits starting and ending at the origin r = form 
families of continuously degenerate orbits, since they can be rotated about arbitrary 
finite angles around the origin due to the U(l) symmetry of the system. This has 
consequences for their semiclassical amplitudes, as discussed in section 14.11 



3.2. Geometry of an arbitrary orbit with v reflections and winding number w 

We now study the geometry of an arbitrary orbit, as illustrated in figure [U Let the 
particle start from a point r (black dot) on a radius vector < r < R (vertical dashed 
line). Let (3 be the starting angle with respect to that radius vector. Let j3' be the final 
angle at which the particle returns to the point r' (circle) on the same radius vector 
after v reflections at the boundary and winding w times around the origin. Both /3 and 
(3' are chosen to lie inside the polygon formed by the orbit. The relation between j3, (3' 
and the reflection angle a at the boundary then is 

= (1 - w)2ir+ (v- l)ir-2va-p. (27) 

Note that according to our definition, all angles a, (3 and j3' are non-negative. Moreover, 
they are restricted in all cases by < a < n/2 and < (3, (3' < ir. The lengths r and r' 
can be expressed in terms of these angles and R by 

sin a , , 

sin [3 

. „ sin a , , 

r' = R—-T,- 29 
sin /r 

§ This is different in the trace formula for the level density: there one integrates over all r and thereby 
includes automatically all pairs of librations starting in opposite directions, so that librations must only 
be counted once 
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Figure 1. Geometry of an arbitrary orbit in the circular billiard with radius R, 
starting at distance r (dot) and ending at distance r' (circle) from the center along the 
same radius vector (shown by the vertical dashed line). 



The condition for an orbit to be closed, r = r', leads to sin (3 = sin j3' which has 
only two solutions in the physically interesting domain j3,(3' G [0, ir], namely: 

a) P' = (3, giving in general a NPO (which, however, can become periodic in an isolated 
point where j3' — (3 — 7r/2), and 

b) j3' — 7r — /3, always giving a PO. 

From (1281 we can express (3 implicitly in terms of r and a by the following relations 



R r 2 1 r 2 

cos f3 — ±— \ — - — sin 2 a , cot (3 = ± \ — - — sin 2 a . (30) 

r V R sin a V R 

The signs above must be taken to be positive if (3 < n/2 and negative if (3 > tt/2. 
3.3. Non-periodic orbits 

Inserting (3' = [3 into ( 1271) . the relation between (3 and a becomes 

P = (1 - w)ir + {v - 1)tt/2 - va . (31) 
Inserting (JHU) into ([28]) yields T clS 8b function of a: 

r(a) = (-l)«i2— ^ r, (32) 

cos(i>7r/2 — va) 

which in general cannot be inverted explicitly. In order to find all NPOs in the circular 
billiard, it is sufficient to look for all real solutions a(r) of the inverse of equation ( |32l) in 
the physical range < a < tt/2. Upon variation of r G [0, R] for a fixed pair (v, w), U(l) 
symmetry breaking at r = creates non-radial NPOs, while bifurcations of some NPOs 
occur at specific points r b > 0, whereby new NPOs or POs are created. It turns out that 
the radial NPOs and the diametric PO are the prime generators of many other closed 
orbits in the circular billiard. The remaining NPOs are created in tangent bifurcations. 
The bifurcation scenarios are discussed in detail in section 13.61 
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The NPOs are symmetric about the radius vector containing the point r. They 
are doubly degenerate due to the time reversal symmetry, which here is identical to the 
reflection symmetry, except for v = 2, w — 1 in the limiting case r = (a = 0, (3 = 7r/2), 
in which case the orbit becomes radial and is mapped onto itself by time reversal. 

The length of the orbit (v, w) becomes 

L V)W (r) = 2(vRcosa + r cos/3) . (33) 

The explicit dependence of (|33|) on v and w is obtained when inserting (13TT) for (3. Note 
that, because of the implicit r dependence of a and (3, the formula ( 133]) does not reveal 
the full (non-linear) dependence on the radius r. 

At r = R, the NPOs become periodic with v + 1 or v — 1 reflections (see 
the systematics in section 13.61) . In this case we have (3 = a or (3 = tt — a, and 
L> v ,w( r = R) = 2R(v ± l)cosa. In the notation of Balian and Bloch [10], the half 
polar angle covered between two successive reflections of an orbit with v reflections is 
4>v,w = wn/v. The lengths then become L vw {r = R) = 2R (v ± 1) sin0„ ±l u ,. 

It is important to realize that the function L VyW (a) given in ( 133]) is stationary with 
respect to a, at a fixed value of r, precisely for the NPOs obeying the relation ( 128]) . 
Indeed, with the extra condition ( T3TT) which implies d/3/da = —v, we find the stationary 
condition to lead to 



dL, 



[a 



da 



2v (r sin (3 — R sin a) = 



(34) 



which is identical with (T28]) . For later reference, we give here also the second derivative: 
d 2 L r 



[a 



=-2ui2cosa(l + v4— )■ (35) 
dcr \ R cos a / 

For the calculation of the Jacobian = J~ V)W defined in (122]) . we make a small 
variation of the initial momentum by a small variation 5(3 of the starting angle. The 
reflection angle then will be a' = a + 5a. From (128]) with fixed r, we obtain 

r cos/5 



5a 



.R cos a 

The particle will then return to the point r + 8r' with the angle /3' 

5(3' = -8(3 - 2v8a . 

From f[2"9"]) we obtain the increment Sr' as 

R 



8(3' , where 



8r' 



sin (3 



(cos a 5a — sin a cot /? 5/3') 



(36) 



(37) 



(38 



Noting that the change of the starting transverse momentum is 8p± = p\8(3 and the 
change of the final transverse distance is 8r\ = Sr' sin (3, we obtain the Jacobian J v w as 



R_ 

Px 



sin a cot (3 + — cos (3 
R 



Using (128]) we can rewrite it in the form 

yJv.w \T 



2v tana cot (3) 



2V R 

— cosp 
Px 



r cos (3 

1 + v- 

R cos a 



(39) 



(40) 
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Like for the lengths (1531) . the explicit dependence of (1591) and (1I0"|) on t> and w comes 
through (I3ip . and their full r dependence is not revealed. At r = R we obtain the values 

2R 

Jv,w(r = R) = — (v ± 1) sin 0„± 1)IU . (41) 

The Jacobian J^, m (r) becomes zero, and therefore the amplitude .A,,,™ defined by 
( 122]) diverges, at three types of critical points: 

1) r = 0: this corresponds to the caustic point of the families of degenerate orbits. 
Making r > leads to U(l) symmetry breaking. 

2a) cos/3 = 0: this corresponds to a point for which the starting and ending angle is 
f3 = 7r/2; the orbit then becomes periodic with v reflections. This point is the 
caustic point — Rcoscj) VtW of the PO family (v,w,l, an isolated member of 
which is created there in a pitchfork bifurcation (see section 13~6]) . 

2b) The term in square brackets is zero, implying 

v r cos (3 = — R cos cu . (42) 
This is fulfilled in the following cases (cf. section [3T6]) : 

a) for the radial type "+" orbits {(3 — ir, a = 0) with (v,w) = (2k + 1, k), k > 0, 
which we call (cf. section [33!) in the points = R/v = R/(2k + 1) at which 
non-radial NPOs (2k + 1, k) are created in a pitchfork bifurcation, or 

b) for non-radial orbits which are born in pairs (v, w), (v, w)' in a tangent bifurcation 
at a critical point r VtW which corresponds to a minimum of the function r(a) in fT28|) . 
This is found from the slope function 

., . dr(a) R r ■ n m /,m 

r (a) = = — 7) — cosasmp + v sin a cos p . (43) 

da sin"' (3 



The condition r'(a = 0) leads, indeed, to (T42l . 

3) r = R: this happens exclusively for the primitive "+" orbit (1,0) which is 
responsible for the Friedel oscillations near the surface [2|. 

The regularizations of the diverging amplitudes, using standard uniform approximations 
for symmetry breaking and bifurcations, will be discussed in section 14.21 

The momentum mismatch needed to obtain the density 8r\(r) in (fT6l) becomes 

Qv,v,{r) = P • P/P 2 = ~ cos(2/3) . (44) 

We stress that the results ( 1331) . ( 1401) and ( 144]) also hold for the radial orbits with 
a = and f3 = or ir, which are discussed specifically in section 13.51 and for which we 
have the relation v — 2k + 1. 
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3.4- Periodic orbits 

For j3' = 7T — /3, we obtain POs. At any point r > 0, they are isolated. If (3' = j3 = tt/2, 
they are symmetric about the radius vector; if (3 ^ 7r/2, they exist in two degenerate 
versions related to each other by reflection at the radius vector. Each of them is doubly 
degenerate due to time reversal symmetry, except for the diameter orbit v — 2. 

Inserting j3' = ir — (3 into fl27|) . the angle j3 falls out and we obtain the reflection 
angle a VyW of the POs: 

( V — 2w\ 71 W7V 

a = a ^ = [-^rr = 2~—> (45) 

and hence 

7T urn 

9v,w = t: - &v, w = — , (46) 

_ V 

as given by Balian and Bloch [10]. The lengths of the orbits which are, of course, 
independent of the starting point r, become 

Ll° w = 2vRsin<p v , w . (47) 



Note that this formula follows from ( T33I) with ( 1461) and cos (3 = 0. 

For the calculation of the Jacobians, we again use equation fl36l) . but the equation 
( |38l) becomes for POs, with sin (3' = sin (3 and cos (3' = — cos (3, 

R 

5r' = (cos a 8a + sin a cot (3 513') , (48) 

sin (3 

where 8(3' is given in fl37|) . With this, we obtain the Jacobians 



rr ( NPO 2V 1 ( 2 d2 • 2 \ 2 ^ (r 2 - i? 2 COS 2 0^,«,) f.tf, 

Jv,w{ r ) = —5 [r-R sm a) = — ^- . 49 

Rpx cos a itpA sm0^ iiu 

The negative sign in front of the last expression above is opposite to that given in |llj . 
where a trace formula for the level density of the circular billiard has been derived. 
It was, however, compensated in [11] by a different book-keeping of the phases in the 
computation of the Maslov index [see our remark after equation ( l57f) below in section 
O] Comparing (j4ll) and (ji9l. we note that J^°(r = R) = -j£$ tW (r = R). At the 
caustic points given by 

= Rcoscj) VtW , (50) 

the Jacobian becomes zero. These are, in fact, the critical points at which all non-radial 
POs (v,w) are created from NPOs (i>,u>) by a pitchfork bifurcation (see section IXBl . 

Trivially, the momentum mismatch is Q v , w {r) = 1 for all POs at all points except 
the reflections points on the boundary. At these, Q V , W (R) can be evaluated as the limit 
r — > R of the expression ( 14*41) for the (suitably chosen) NPOs. 
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3. 5. Radial orbits and bifurcations of POs 

There are only two kinds of radial NPOs starting at any point < r < R in the radial 
direction. They correspond to the only NPOs that exist in one-dimensional systems. 
Like in PQ, we call the "+" orbits those which start in the outwards direction and 
are reflected at the nearest turning point, and the "— " orbits those which start in the 
inwards direction and are reflected at the opposite turning point. The primitive "+" 
and "— " orbits have only v = 1 reflection at the boundary; by definition their repetition 
number is k = 0. Using the geometry of section [3T3l we have a = and (3 = ir or (3 = 
for the "+" or "— " orbit, respectively. To both orbits, one may add k > full librations 
between the two opposite turning points; we call k their "repetition number". Their 
reflection number then is v — 2k+ 1. In order to be consistent with fl3TT) . we define their 
winding numbers to be w = k for the "+" and w = k + 1 for "— " orbits. Because both 
v and w are determined by k, the repetition number k is sufficient to characterize the 
radial NPOs uniquely. We shall denote them here by the symbol L± . At r = R, the 
orbits with k > become identical with the fc-th repetitions of the primitive PO 
(2,1), while the orbits L_ for all k become its (k + l)-th repetitions (2k + 2, k + 1). 

The formulae ( 1331) . ( )40l) and ( )44l) can be used with cos/3 = =Fl for the L± orbits. 
One of their characteristic features is that they have opposite initial and final momenta, 
p' = — p, so that Q±\r) = — 1 at all points. Their lengths become 

Lf(r) =2[(2k + l)RTr}. (51) 

For the Jacobians we obtain 

ji k \r) = ^-r[(2k+l)r T R]- (52) 

Rpx 

Both Jacobians are zero at r = where the "+" and "— " orbits become identical 
and form a rotationally degenerate family, as discussed in more detail in section 14.11 
Furthermore, the L+ orbits have zero Jacobians at = R/{2k + 1). For the primitive 
"+" orbit, this gives the turning point r = R at the boundary, where the divergence 
of the semiclassical amplitude is removed by the uniform approximation discussed in 
section 14731 At the points rj, = R/(2k + 1) with k > 0, non-radial NPOs with 2k + 1 
reflections at the boundary are born from the orbits in pitchfork bifurcations. The 
example of k — 1, where the new NPO here is called the A orbit (see also figure [2]), will be 
discussed explicitly in section 14.21 along with the corresponding uniform approximation 
for its semiclassical amplitude. 

3. 6. Creation and systematics of non-radial orbits 

The properties of all closed orbits in the circular disk are obtained by systematically 
inverting fl32|) for all (v,w) and inserting the resulting angles a(r), /3(r) at each point 
r into the equations fl33l) . fT40l) and fj44"|) . As already mentioned, bifurcations of the 
orbits occur under variation of r, at which new orbits are created. Furthermore, the 
U(l) symmetry of the continuously degenerate orbit families existing at r = is broken 
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when r becomes > 0, whereby non-radial NPOs are created. In fact, we find that all 
non-radial NPOs are created either (i) by the k-th repetition of the diametrical PO 
(v, w) = (2k, k) (with k = 1, 2 . . .) by U(l) symmetry breaking, or (ii) by the orbits 
by pitchfork bifurcations, or (Hi) pairwise by tangent bifurcations. All non-radial POs 
are created by pitchfork bifurcations from non-radial NPOs. 

The detailed systematics are as follows. (We include here also the Morse indices 
whose determination is discussed in section I3TT1 ) 

• Non-radial orbits have v > 2 and 1 < w < |_| J , where |_^J is the integer part of x. 

• NPOs with even v and maximum w, i.e. (v, w) = (2k, k), are created from the family 
of POs (2k, k) (i.e., the k-th. repetitions of the diameter orbit) by U(l) symmetry 
breaking at r = which is the caustic point of the family. These new NPOs exist 
at all distances < r < R. Their Morse index is /i = 3v . At r = R they become 
equal to the POs (2k + l,k). Analytical expressions for the orbit (2,1), which we 
here denote by V, are given in Appendix A. 1[ 



• NPOs with odd v and maximum w, i.e. (v, w) = (2k + 1, k), are created from the 
radial orbits by a pitchfork bifurcation at r k = R/(2k + 1) which are critical 
points of (r). They exist at all distances r k < r < R; at r = R they become 
equal to the POs (2k + 2,k). At r = they create the PO (2k + 1, k) in a further 
pitchfork bifurcation, whereby is the caustic radius ( 1501) of that PO family. The 
Morse index is fi = 6k + 2 for r k < r < and \i = 6n + 3 for < r < R. 

• All remaining NPOs are created in pairs (v,w), (v,w)' (with identical values of v 
and w) by tangent bifurcations at critical points r VjW (with < r VtW < R) where the 
Jacobian J VyW becomes zero [cf. equations ( H3l) . (H2"l) and (HQ])]. They exist only for 
r v,w < r < R. The points r v<w occur as minima of the function r(a) in fl32l) . The 
two orbits correspond to the two branches of the inverse function a(r), whereby the 
orbit (v,w) has the slope r'(a) > and the orbit (v,w)' has the slope r'(a) < 0. 
They have the following properties and behaviors: 

- The orbit (v, w) creates the PO (v , w) by pitchfork bifurcation at the caustic 
radius > r v,w Its Morse index is [X — 3v — 1 for r v , w < r < and /i = 3v for 
r v!w < r < R. At r = R, the NPO (v, w) becomes equal to the PO (v + l,w). 

- The orbit (v, w)' undergoes no further bifurcation. Its Morse index is /x = 3v — 2 
at all radii r V)W < r < R. At r = R, the NPO (v, w)' becomes equal to the PO 
(v — l,w). 

In figure [2] we illustrate some of the bifurcation scenarios. In the upper panels we 
show the lengths L v>w (r) and in the lower panels the Jacobians J V)W (r) of some orbits as 
functions of the starting point r. In the left panels, the radial orbit is seen to undergo 
a pitchfork bifurcation at r = R/3, at which the non-radial orbit A (v,w) = (3,1) is 
created. At r = R/2, the orbit A undergoes a pitchfork bifurcation at which the 
triangular PO (3, 1) (noted as A) is born. In the right panels, we see the creation of 
a pair of orbits P (4,1) and P' (4,1)' by a tangent bifurcation at r 41 = 0.682489-R, to 
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Figure 2. Bifurcations of the NPOs (u, w) = (3,1) (L4/ and A) and the triangular 
PO (3,1) (A) (left panels), and of the NPOs (v,w) = (4,1) (P), (4,1)' (P') and the 
square PO (4,1) (□) (right panels). Top: lengths L VyW (r), bottom: Jacobians J v , w {r). 
In the lower panels, the subscripts of the orbit symbols indicate their Morse indices. 






Figure 3. Shapes of the NPO A (3,1), shown by the solid lines, at the starting points 
r = 0.345-R, r = R/3 and r = Q.75R (from left to right). In the right panel, the 
dashed lines show the shape of the triangle PO A, which exists in two versions lying 
symmetric to the diameter that contains the starting point r (black dot). 



the left of which they do not exist. The orbit P undergoes a pitchfork bifurcation at 
= R/y/2, at which the squared PO (4,1) (noted as □) is born. In the lower panels, 
the subscripts of the orbit symbols indicate their Morse indices \x. 

In figure [3] we show the shape of the A orbit at three starting points (from left 
to right): at r = 0.345i? shortly after its bifurcation from the orbit, at r = R/3 
(i.e. at the pitchfork bifurcation where it is identical with the PO A), and at r = 0.75i? 
together with the two POs A passing through the same point. Three shapes of the orbits 
P and P' are shown in figure HI chosen (from left to right) at the starting point r = r^i 
(tangent bifurcation, where they are identical), at r = R/y/2 (pitchfork bifurcation, 
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Figure 4. Shapes of NPOs P (4,1) (solid lines) and P' (4.1)' (dashed lines) with 
starting points r = n t i, r = R/V2 and r = 0.8R (from left to right). 



where P is identical with □), and at r = 0.8-R. Analytical expressions for the lengths 



and Jacobians of the orbits A and P, P' are given in Appendix A. 2 and Appendix A. 3 
respectively. 




Figure 5. Small perturbations (thin dashed lines) of an orbit (thick solid line) starting 
at a point r (black dot) with distance sq to the first reflection point. The first conjugate 
point is shown by a circle. Left: sq > 2-P/3; here the conjugate point occurs at distance 
si from the first reflection point. Right: sq < 2P/3; here the conjugate point occurs 
only after the second reflection at distance t\ from the second reflection point. 



3. 7. Calculation of Morse indices 

Following Gutzwiller [3], we calculate the Morse index of a given orbit by counting 
the number of conjugate points along the orbit. These are by definition the points r c 
in which a fan of infinitesimally perturbed orbits, all starting at the same point Tq, 
intersects itself. To this number, one has to add twice the number of reflections at the 
boundary, since each hard-wall reflection yields a phase shift of it, while each conjugate 
point that does not lie on the boundary yields a semiclassical phase tt/2 [3] (except 
for multiple conjugate points which can occur in the presence of higher continuous 
symmetries in D > 2; this does not happen in the present system). 

We illustrate the procedure of determining the conjugate points in the circular 
billiard in figure [5] for an arbitrary non-radial orbit; the resulting formulae ( 1551) . ( |56l ) 
hold, however, also for the radial NPOs. Let the orbit start at a distance r from the 
center. Let the starting angle be (3, as in figured], and denote the distance of the starting 
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point to the first reflection point by Sq. Let the distance of two successive reflections of 
the orbit be 2P, so that we have 

P = R cos a . (53) 

We now consider a fan of orbits, starting at the same point r, with infinitesimally 
perturbed angles (3' = (3 + 5(3. These perturbed orbits will intersect in a conjugate 
point, shown in figure [5] by a circle. For sq > 2P/3, the first conjugate point is located 
on the next portion of the trajectory at a distance si from the first reflection point, 
as shown in the left panel of figure [5j For so < 2P/3, the first conjugate point does 
not lie on the next portion of the unperturbed trajectory, but on the subsequent one, 
i.e., between the second and third reflections, at a distance t\ from the second reflection 
point (see the right panel of figure [5]). 

The calculation of the distances s\ and ti is elementary, using differential calculus 
and trigonometry. The perturbed reflection angle at the boundary is a' = a + 5a with 

. r cos(3 . . 

Sa=- -5(3, (54) 

R cos a 

as follows from fl28|) and (|3~T1) . For infinitesimal perturbations {5(3 — * 0), s\ and t\ must 
be independent of 5(3 due to the definition of a conjugate point. This leads, after some 
algebra using (153]) and r cos (3 = Sq — P, to the following values which are consistent 
with the findings of |12j : 

3 ^ p -j^i for s » a I p - (55) 

Special cases are the following: (i) s = 2P/3 =>- si = 2P or t\ = 0. (ii) s = 
=> h = 2P/3. (Hi) s = 2P Sl = 2P/3. (iv) s = P => s x = P. In this last 
case, the starting point is a caustic point of a family of degenerate orbits due to the 
rotational U(l) symmetry. We note that these caustic points, common to families of 
rotationally degenerate (periodic or non-periodic) orbits, never coincide with conjugate 
points unless they are chosen as starting points (which, however, one should avoid as 
emphasized below). 

The successive conjugate points s n or t n with n = 2, 3, . . . along an orbit are 
then found iteratively by taking the previous conjugate points as new starting points, 
evaluating their distance s' Q to the next reflection point, and using the relevant formula 
of ( 1551) . (l56l) with so replaced by s' . 

As stated above, the total Morse index \x for a given orbit is obtained by adding to 
the number of conjugate points twice the number v of reflections at the boundary. We 
emphasize that for the resulting Morse index to be well defined, the starting point may 
not he any critical (caustic, conjugate or reflection) point (cf. section H]). 

For the POs (v,w), we find in this way that the number of conjugate points is 
always v — 1, i.e., one less than the number of reflections. Including the contributions 
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Figure 6. Orbit L^ 1 ' with starting point at r = 0.15R (left) and r = 0.75R (right). 
Fans of small perturbations with different starting angles exhibit the positions of the 
conjugate points. The total Morse indices become // = 6 + 2 = 8 (left, with two 
conjugate points) and fx = 6 + 1 = 7 (right, with one conjugate point), in agreement 
with the general result (158)) . 




Figure 7. Same as figure [6] for the orbit A with starting point at r = Q.35R, giving 
/j = 6 + 2 = 8 (left, with two conjugate points) and r = 0.75i? giving ^ = 6 + 3 = 9 
(right, with three conjugate points). 



from the reflections, their total Morse index therefore is given by 

H™ = 3v-1. (57) 

This Morse index has been used implicitly in the derivation of the trace formula for the 
level density of the circular billiard in [UJ. The calculation there was, however, done 
with a different book-keeping of the phases, absorbing the negative overall sign of J v<w 
in our result (1491) into the Maslov index given in [TT] . 

For the radial NPOs we obtain the following Morse indices: 

$ = 6^ + 2 for r <r k = R/(2k + 1) , 
fjM = 6k + l for r>r k = R/(2k + 1) , 

fjW = Qk + 3. (58) 

For the non-radial NPOs, no single general formulae can be found. We have, however, 
given the Morse indices of the various types of NPOs in their systematics at the end of 
section 13.61 

The geometric construction of the conjugate points for the orbit is illustrated in 
figure O and for the orbit A, whose bifurcation is shown in figure |2] and further discussed 
in section 14.21 in figure [71 
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3.8. Relation to multiple-reflection expansion of Green function 

In a publication that appears to have remained largely unnoticed by the semiclassics 
community!]]] Bonche [13] has applied the multi-reflection expansion of the single-particle 
Green function developed by Balian and Bloch [10] for calculating particle densities 
of three-dimensional billiards. In particular, he gave explicit results for the spherical 
billiard whose closed orbits are identical to those in the two-dimensional circular billiard, 
although their semiclassical amplitudes differ since they depend on the dimension D. 
Bonche discussed the orbits with one reflection, i.e., the orbits in our notation, and 
obtained similar results to ours for their regularized contributions at the center and the 
boundary (cf. sections 14.1.11 and 14.31 below). For the orbit with two reflections, i.e., the 
orbit V (2,1) in our notation (cf. Appendix A.l ), he obtained the result (IA.2I) but gave 



an obviously wrong expression for its length ( 1A.3I) that would become zero for r = R. 



Bifurcations of NPOs have not been noticed in [T3l. 



| We are grateful to A G Magner for drawing our attention to this paper which he discovered recently 
in the preprint collection of the late V M Strutinsky. 



Closed orbits and spatial density oscillations in the circular billiard 



19 



4. Regularizations of semiclassical amplitudes 

At the critical points of ^J v ,w{f) discussed in the previous section, the semiclassical 
amplitudes Ay (A, r) in the equations f|T4l) - ffTB"]) for the spatial densities diverge. The 
reason for these divergences can be traced back to a break down of approximate stat- 
ionary-phase integrations used in the derivation of Gutzwiller's semiclassical Green 
function [3] that underlies also our semiclassical theory of spatial density oscillations [I]. 
In order to regularize the amplitudes, we resort to standard uniform approximations that 
have been developed for semiclassical trace formulae for the level density in terms of POs 
in connection with symmetry breaking [HJ [15], [16j [17] and bifurcations [T4" l [TB I [T8 | [19 ] 120] . 

For symmetry breaking, the uniform approximations in general depend not only 
on the particular symmetry that is broken, but also on the way in which a system 
is perturbed [T71 I2T] . For U(l) symmetry breaking, however, the uniform approxi- 
mation developed in [151 EH] is universal and can be readily applied in section 14.11 
For bifurcations, on the other hand, the uniform approximations do not depend on 
the particular system. They are, in fact, universal for each generic type of bifurcation, 
making use of the standard normal forms known from bifurcation theory (see, e.g., |22j). 
It turns out that the pitchfork and tangent bifurcations occurring for the NPOs in the 
circular billiard have exactly the same generic features as those of POs, so that we readily 
can make use of the corresponding uniform approximations which were developed in [18J. 
That for pitchfork bifurcations will be discussed in section I4~2l Finally, the contribution 
of the primitive "+" orbit at r = R to the density p(r) can be regularized as already 
stated in [2] and briefly illustrated in section 14.31 All these regularizations will then be 
used in our numerical calculations in section [5j 

We shall only present explicitly the uniform approximations for the density 
oscillations Sp(r), whereby the basic input is given by the amplitudes A 7 (r) defined in 
( 1221 . which diverge at the critical points. For symmetry breaking and for bifurcations, 
the uniform approximations for the kinetic energy densities are then easily obtained 
according to the semiclassical equations ( fTBT) and ( fT6l) : For Sr(r), we just have to 
multiply the result of Sp(r) by the factor p\/2m. For 7"i(r), we must furthermore replace 
the amplitudes A 7 (r) in the resulting expressions for 5p{r) by ^4. 7 (r)Q 7 (r). This does, 
however, not work for the regularization at the boundary r = R, as discussed at the end 
of section 14.31 

4-1- Symmetry breaking at r = 

4- 1.1. Radial NPOs L±^ . The local regularization for the radial orbits at r = has 
already been given in [TJ Ej and yields the result (1251) which is valid for r <C R. Here we 
first present an alternative derivation of this result, following the perturbative approach 
of Creagh [TTJ, using r/R as the perturbation parameter which for r > breaks the 
U(l) symmetry. We then follow the approach of [T5l [16] to derive a global uniform 
approximation that is valid up to all distances r for which no new critical phenomenon 
occurs. (We refer to section 6.3.1 of [5] for an easy understanding of the main ideas and 
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the technical details.) 

The basic idea is that the NPOs L± at r > can be considered as the result of a 
symmetry-breaking process. At r = 0, they are degenerate, forming a family of orbits 
with U(l) symmetry obtained by rotation about an angle £ [0,2tt). For r > this 
symmetry is broken; the scenario corresponds to the generic case of breaking a U(l) 
torus into a pair of isolated orbits according to the Poincare-Birkhoff theorem [22J . 

We make use of the stationary property of the length function L vw (r) in fl33l) at a 
fixed distance r and for fixed values v,w. With the help of ( 128]) we rewrite it as 

L V)W (r, j3) = 2v\J R 2 - r 2 sin 2 (3 + 2r cos (3 . (59) 

At the stationary points j3_ = and j3 + = 7r, this yields exactly the lengths of the orbits 
L± given in fl5"Tj) with v — 2k + 1. For r = 0, the function (159]) is independent of (3 and 
yields the lengths 2(2k + 1)R of the degenerate orbit families, in short called "tori". 

We now re-interpret the L v w (r, (3) in the following way. We rename the variable 
f3 to and take it as the parameter £ [0, 27r] describing the U(l) symmetry of the 
tori. For r > 0, L VtW (r, 0) yields the lengths of the perturbed tori. To first order in the 
perturbation parameter r/R, we thus get the actions of the perturbed tori to be 

S {k \r,<j)) = S { h) + 2rp A cos0, S { Q k) = 2p x {2k + 1)R . (60) 

According to the classical perturbation theory of Creagh [17], we can now write the 
contribution of the orbits to the particle density in ( fl4l) as the integral 

StP ( r) = Re ±Y / i fc ( r ,0) e ^^ ) M)M-^ ) V2-3./4] d0; (61) 

(k) 

where fi are the Morse indices of the unperturbed tori which will be determined below. 
The amplitude functions Ak(r, 0) will also be determined in the following; they must be 
chosen such that for r ~ the sum of integrals ( 16TI) leads to the result ( 1251) . 

By construction, the stationary points of S^ k \r, 0) at fixed r are 0_ = and + = n 
and yield the actions of the isolated Lj9 orbits: 

S {h \r, 0±) = S { ±\r) = 2p x [(2k + 1)R t r] . (62) 

Hereby + and 0_ are the directions into which the orbits start from the point r > 0: 
the "+" orbits start towards 0+ = it, while the "— " orbits start to the opposite direction 
0_ = 0. The stationary-phase approximation for the integrals ( 16TT) . with appropriately 
chosen values of v4^(r, 0), therefore yields exactly the contributions of the isolated 
orbits to (TT^j) . Since S"( k \r,(j) + ) > and S"( k \r,(j)_) < 0, their Morse indices differ by 
one unit as they must, according to their values given in (158]) for small r. 

We now take, in a first step, the amplitudes A k (r, 0) to be independent of 0, using 
their forms (122]) valid for small r, and perform the integral in flBTj) exactly rather than by 
stationary-phase approximation. This leads to the Bessel function Jo(2rpx/h) with the 
correct argument, hence reproducing the locally regularized density ( 1251) after summing 
over all repetitions. This result corresponds to a local uniform approximation in the spirit 
of [TU [T7] which is correct locally for r ~ for which the amplitudes were approximated. 
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(We recall that using the full amplitudes given by (1221) . the summation over k cannot 
be done analytically.) 

In the second step, required to obtain a global uniform approximation valid also 
for large r, we choose A k (r, <fi) to be of the form A k (r,(p) = a k {r) + b k (r) cos<p. The 
coefficients a k {r) and bk(r) are chosen such that they reproduce the exact amplitudes 
A± (r) given in fl22|) at large distances r from the center, where the stationary-phase 
approximation for the integral is valid. The exact integral for the terms including b k (r) 
leads to the Bessel function Ji(2rp\). As a result, we obtain the following global uniform 
approximation for the combined contributions of the radial NPOs L± to the particle 
density oscillations: 

OO s 
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Hereby (r) are the exact amplitudes obtained from (1221) (we have omitted the Fermi 
energy A in their arguments), using the properties of the orbits given in section [3751 
The result f l6"3"|) is finite at r = 0, where the second term oc J\(2rp\) does not 
contribute. This term becomes important, however, at larger distances r. When 
2rp\ ^> h, we may use the asymptotic forms of the Bessel functions for large arguments 
[23] (which is tantamount to using the stationary-phase approximation for the above 
integrals), and then the result ( l63l) goes over into the contributions of the isolated 
orbits corresponding to ([T 
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In order to obtain the correct Morse indices ( 1581) (valid below the bifurcation points r k of 
the "+" orbits), we now see that we have to choose = 6/c + 5/2, which corresponds to 
the average value of and [/,_ . With this, the phases can be simplified to a common 
sign factor and the result (l63l becomes 



5 r p 



(un). 



V^rp x /h J2(- 1 ) k {Mr)Jo(2rp x /h) cos[S^/h] 



k=0 



-AA k (r) Jx{2r Px /h) sin[Sf/fr]}. (66) 

The result f l66|l can only be used below the critical points r k = R/ (2k + 1) at which the 
"+" orbits with k > 1 bifurcate, or below the boundary r = R where the amplitude 
for the primitive "+" orbit with k = diverges. Near these critical points, we need 
different uniform approximations that will be discussed in sections 14.21 and 14.31 
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Figure 8. Non-radial NPO (2,1) (triangle, in the text called V) starting from a point r 
(black dot) near the center, where it was created from the PO (2,1) (vertical diameter) 
by U(l) symmetry breaking. The arrows indicate the starting directions of the four 
isolated orbits (see text for details). 

4.1.2. Non-radial NPOs (2k, k). The non-radial NPOs with (v,w) = {2k, k) and 
k = 1,2,... are created from the k-th repetitions (2k, k) of the diametrical PO by 
symmetry breaking at r = (see the systematics in section [3T6T) . For small r, they leave 
their starting points r perpendicular to the diameter containing the POs (2k, k). As 
an example, we show in figure M the NPO (2,1), which we here call V, with a starting 
point r (black dot) close to the center. For r = 0, these NPOs are degenerate with 
the families of diametrical POs and their semiclassical amplitudes diverge. The uniform 
approximation for these NPOs can be derived in a similar way as above for the radial 
NPOs. We first give the derivation for the orbit V (whose analytical properties are 
given in the appendix) and then generalize the result to all NPOs (2k, k). For r = 0, the 
PO orbit family (2, 1) has the action Si T = &p\R- For r > 0, the family breaks up into 
a degenerate pair of isolated diametrical POs, starting in opposite directions along the 
diameter going through r (this is the vertical diameter in figure [8]), plus a discrete pair 
of NPOs V related by time-reversal symmetry (cf. the triangle in the figure, starting in 
the two directions indicated by the arrows). Their action is S 2 (r) = p\L 2 ^(r) with the 
length .£,2,1 (r) given in (1A.3I) . We can obtain the actions of all four isolated orbits from 
the stationary points of the action function 

S 2 (r, 0) = S lr + 2AS(r) sin 2 , AS(r) = r 2 p x /R . (67) 

At (j) = and ir, it yields the action S\ r of the diameter orbit, and at <fi± = ±7r/2, it 
yields the perturbed action 5i r + 2AS(r) of the V orbit, valid for small r as given in 
(1A.5j) . Note that the values of <fi at these stationary points correspond to the directions 
(shown by the arrows in figure EJ) in which the perturbed isolated orbits leave the point 
r (for the V orbits only in the limit r — > where the symmetry breaking occurs). 
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We thus write the common uniform contribution of the four perturbed isolated 
orbits to the particle density (JUj) as the following integral: 

1 



Sp 2 ,i(r) = Re 



2ty 



A{r, 0) e i[&(r^)/a- M7 r/2- 37 r/4] ^ _ 
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For a constant amplitude A(r, 0) independent of 0, the integral in ( 1681) yields the Bessel 
function Jo(AS(r)/h), corresponding to a local uniform approximation valid near r = 0. 
Adding to A(r, 0) a term proportional to sin 2 yields an additional Bessel function J\ 
with the same argument. Proceeding like above for the radial NPOs and replacing the 
perturbative action change AS(r) in the arguments of the Bessel functions by the exact 
difference AS(r) given in (J7TT) below, we find the global uniform approximation for the 
combined NPOs and POs (2,1) to be 
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The amplitudes above are given by 
Mr) = l[Au(r)+A 2 ,i(r)] 



AA(r) 



[4(r)-AiW 



(69) 



(70) 



where «4i r (r) is the amplitude ( 1221) evaluated for the diameter PO (2,1) in term of its 
properties given in section [33 and «4 2 ,i ( r ) that of the V orbit using its exact properties 
given in (1A.3I) . (1A.4I) . Both amplitudes must include a discrete degeneracy factor 2 
corresponding to their two starting directions. Furthermore, the quantities AS'(r) and 
S{r) are defined by 



AS{r) = - [S 2>1 {r) - Sy , S{r) = - [S 2>1 (r) + S lr ] 



(71) 



where S 27 i(r) is the exact action of the V orbit, given by 5*2,1 (r) = P\L2,i(r) using ( 1A.3I) . 

For large enough distances r such that AS(r) ^> h, we can use the asymptotic 
forms of the Bessel functions for large arguments, and (169!) goes over into the sum of 
the contributions of the isolated diameter POs and V orbits 
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(72) 

5 for the diameter 



With the choice /io = 11/2, we obtain the correct Morse indices /ii r 
orbit and /i 2; i = 6 for the V orbit. 

The uniform result (1691) gives a finite combined contribution of the degenerate orbit 
families at r = 0. Since the orbit V undergoes no bifurcation and its amplitude stays 
finite for all r > 0, the uniform result ( |69l) is valid in the entire region < r < R. 

The result ( |69l) looks very similar to the uniform contribution (166|) of the radial 
NPOs, except that the argument of the Bessel functions here is quadratic in r (for small 
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r), which yields a slower departure from its maximum. This is characteristic of the 
long-ranged irregular oscillations Si TT p(r) that are due to non-radial NPOs. In fact, the 
presently discussed class of non-radial NPOs (2k, k) is dominating the irregular part of 
the density oscillations near r = 0, as we shall see in section [5j 

We therefore now generalize ( 1691) to include all NPOs of type (2k, k) that are created 
by symmetry breaking at r = from the /c-fold repetitions of the periodic diameter orbit 
(2,1). For that we simply have to substitute the amplitudes ^.2,1 ( r ) and actions S21 (r) 
in the equations ( 17D1) . ( JTTT) by those obtained from the general formulae ( 1331) and ( 1401 ) 
for the NPOs with (v, w) = (2k, k), and the quantities Ai T (r) and Si T by those obtained 
from (T4"Tj) and (H9l) for the POs. The Morse index in ( 1691) must be substituted by 
p^ = Qk — 1/2. After these substitutions we must, of course, sum over all k > 0. 
Taking out a common phase factor, we obtain the global uniform approximation for the 
dominating contribution to the irregular density oscillations: 
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where the amplitudes are 
1 



Mr) = o [4W + A 2k , k (r)] , AA k (r) = - [A kv (r) - A 2k , k (r)] 



and the action functions are 
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The subscripts "/cr" refer to the diametrical POs with their properties given in section 
13. 4[ and the subscripts u 2k, fc" to the NPOs with their properties given in section 13. 3[ 

The NPOs (2k, k) exist at all r > and undergo no bifurcations. The uniform 
result ( 1731) is therefore valid in the whole range r e [0, R}. In the numerical calculations 
presented in section [5J it will be seen that the summation over k can be practically 
limited to k max ~ 10 — 20 in all cases. 



4-2. Uniform approximations for bifurcations 

Here we discuss the uniform approximations needed to regularize the divergences of the 
semiclassical amplitudes arising at bifurcations of the closed orbits. To our knowledge, 
there exists no systematic study and classification of bifurcations of non-periodic closed 
orbits in the literature. We therefore had to convince ourselves that the bifurcations 
occurring in our present system are of the same generic types as those known j22l 124"] 
for periodic orbits. This appears, indeed, to be the case. For the pitchfork bifurcations, 
we shall demonstrate this in a particular example in section 14.2.11 and then present the 
corresponding uniform approximation for the semiclassical particle density. In section 
14.2.21 we briefly discuss the tangent bifurcations. 
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4-2.1. Pitchfork bifurcations. In the left panels of figure [2] we have shown two successive 
pitchfork bifurcations. We shall take the example at the critical radius r = R/3 where 
the non-radial orbit A (3,1) is created from the radial orbit L+ . We draw particular 
attention to the bifurcation diagram on the lower left of figure [2j For isolated periodic 
orbits of D = 2 dimensional systems, the corresponding diagrams are obtained when 
the trace of the stability matrix, TrM(e) as function of a bifurcation parameter e, is 
considered. We assume a bifurcation to occur at the point e = 0. In the trace formula 
for the level density of such systems [3], the quantity Tr M(e) —2 appears under a root in 
the denominator of the semiclassical amplitudes, which diverge at the bifurcation point 
where TrM(O) = 2. Now, it is a characteristic feature [TBI 120] of a generic pitchfork 
bifurcation that the slopes of the functions TrM(e) of the two orbits involved in the 
bifurcation fulfill the following relation at the bifurcation point: 

TrM' B (0) = — 2TrM^(0). (76) 

Hereby A is the parent orbit which changes its stability at the bifurcation point and B 
the daughter orbit that is created at the bifurcation and only exists on one side of it. 
(For a mathematical proof of this "slope theorem", we refer to |25j.) 

In our present system, the decisive quantity under a root in the denominator of the 
semiclassical amplitudes (1221) is the Jacobian j7„ jU ,(r), and r is the bifurcation parameter. 
In the example of figure [2] under discussion, L+ is the parent orbit A and A is the 
daughter orbit B existing only above the bifurcation point r = R/3. Using the explicit 
results ( 152|) and (lA.lOp for the Jacobians of these two orbits, we find 

4 



-R/3 PA dr 



= ~ , (77) 

r/3 Pa 



so that the slope theorem ( 1761) is, indeed, satisfied. 

A pitchfork bifurcation which is generic in the sense of Meyer [23] is period-doubling: 
the new orbit B has twice the period of the parent orbit A at the bifurcation point. 
However, in systems with discrete symmetries the pitchfork bifurcations are typically 
isochronous [201 [2S1 [27]: both orbits have the same periods at the bifurcation point, 
but the daughter orbit B has a double discrete degeneracy compared to the parent orbit 
A. For the NPOs discussed here, the pitchfork bifurcation is of the second, non-generic 
type: their running times T(r, A) (given here simply by their lengths) are identical at 
the bifurcation point, but the A orbit has a degeneracy of two due its time-reversal 
symmetry while the orbit is not degenerate. 

Finally, the following rule is known [18] for the Maslov index appearing in the trace 
formula: at the pitchfork bifurcation, the index of the parent orbit A changes by one 
unit, and the index of the orbit B existing after the bifurcation is the same as that of 
the orbit A before the bifurcation. Using our results in section 13.51 and the appendix, 
we see that this rule is, indeed, also fulfilled here by the Morse indices: the Morse index 
of the orbit changes from 8 to 7 at r = -R/3, and the A orbit has the Morse index 
8 close to r > R/3. After these affirmations of the nature of the present pitchfork 
bifurcation, we can immediately apply the uniform approximation that was developed 
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in [18] to obtain the following common contribution of the and A orbits to the 
particle density oscillations in terms of Bessel functions of non-integer orders!^ 
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(78) 



in terms of the exact amplitudes ^4 3j i(r), A 



and orbits, respectively. The choice of fio = 7 yields the correct Morse indices of 
the isolated orbits on either side of the bifurcation. Note that for r < R/3, there exists 
an imaginary "ghost" orbit A which contributes with real action S^i (r) and amplitude 
-4.3,1 ( r ) which are just the analytical continuations of these properties of the real A 
orbit, as given in (1A.8j) and (lA.lOp . Note also that the amplitude A^\ must include the 
degeneracy factor 2 due to the two time orientations of the A orbit. 
At the bifurcation point, the result (1781) takes the finite value 
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The uniform approximation (1781) is again global in the sense that it goes over into the 
sum of the isolated contributions of both orbits (including the degeneracy of the A orbit) 
sufficiently far away from the bifurcation, where the argument of the Bessel functions is 
large enough so that their asymptotic forms can be used. 

Having confirmed that all other pitchfork bifurcations occurring in the circular 
billiard are of the same type, we can apply the result (178"|) to them, replacing in (1791) 
the amplitudes and actions of the orbit by those of the parent orbit, those of the 
A orbit by those of the daughter orbit (irrespective of whether this is a NPO or a PO 
like in the other pitchfork bifurcations seen in figure [2]), and cr by sign(r — r^) using 
the appropriate bifurcation radius r&. The Morse index /xq in (178]) has to be chosen 
according the rules given in HE 



% We can also derive this result by making an ansatz similar to ((68)) and expanding the action function 
^3,i{r, 0) = p\L3 t i(r, /3) using (|A.8[) in the phase of the integrand around the stationary point r — R/3, 
/3 = ir. Expanding it up to first order in e = r — R/3 and fourth order in 5 = (3 — n, we obtain exactly 
the (one-dimensional) normal form of the action for the pitchfork bifurcation that was used in [18j to 
derive the uniform approximation 
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We must emphasize that the uniform approximation (I78p only holds up to distances 
from the bifurcation point r, until which the respective orbits do not change their 
nature by a new bifurcation (or by symmetry restoring when going towards r = 0). 
This becomes a practical problem when successive bifurcations lie close, i.e., when the 
relevant action difference between the two points is of the order of, or smaller than h. 
We discuss this problem again in the following section and in section [51 

4-2.2. Tangent bifurcations. The tangent bifurcations of POs in systems with D = 2 
always occur as generic isochronous bifurcations pH [26] . The two orbits created thereby 
have infinite slopes of their stability traces [T8l [25] : 

lim TrM'(e) = ±oo . (81) 

The same holds for the Jacobians of the NPOs created by tangent bifurcations in our 
present system, as seen for the example of the orbits P and P' on the right side of figure 
[2j Characteristically, too, the Morse indices of the two orbits differ by one unit close to 
the bifurcation (i.e., before any further bifurcation occurs). 

The global uniform approximation for the tangent bifurcation has been derived in 
[T8] ; it looks similar to ( 178]) except that the Bessel functions have the orders ±1/3 and 
±4/3 (and can also be expressed in terms of Airy functions, cf. also [20]; on that side 
of the bifurcation where the orbits are not real, modified Bessel functions occur). 

The reason for not quoting the result of [18] here is that for all NPOs in the 
circular billiard created in this way, the tangent bifurcation occurring at the critical 
point r = r VjW is immediately followed by a pitchfork bifurcation of one of the orbits 
at a nearby critical point r = > r = r VjW where the PO (v, w) is created (see the 
systematics in section 1376]) . The two bifurcations lie close in the sense that the action 
difference S ViW (r = r^°) ~ S V)W (r = r v>w ) is not much larger than h~. Consequently, the 
arguments of the Bessel functions in the uniform approximations are never large enough 
to allow for their asymptotic expansion. In other words: in the region between the 
two bifurcations, the two orbits are never sufficiently isolated for the two bifurcations 
to be treated separately. This necessitates to treat what in the relevant literature 
is called a "codimension-two" bifurcation. Such bifurcations have been investigated 
and semiclassically applied in [281 EE]. We have, however, not been able to find an 
appropriate uniform approximation for our present bifurcation scenario which is similar 
to, but not identical with the "butterfly catastrophe" discussed in [28]. Luckily, the 
NPOs created by tangent bifurcations are practically negligible (cf. end of section [571]) . 
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Figure 9. Contribution of primitive orbit to the particle density in the circular 
billiard with N = 606 particles {h 2 /2m — R = 1). Solid line: uniform approximation 
(|82[) . dashed line: asymptotic expression (I83[) . 



4-3. Regularization at the boundary r = R 

We have already mentioned in section 13.31 that the semiclassical amplitude of the orbit 
which is responsible for the Friedel oscillations [2], diverges at the boundary 
r = R due to the zero of its Jacobian ( 152|) . A method to regularize this divergence 
for the density using a uniform approximation of the Green function for short times was 
proposed in [30] . The explicit result for D- dimensional spherical billiards has been given 
in [2] ; the particular result for D = 3 was also obtained in [13] . For the two-dimensional 
circular billiard we obtain the following uniform contribution of the primitive "+" orbit 
L<p to the particle densityQ 

2nh V r (R — r) 

Sufficiently far from the boundary, this result goes over into the correct expression for 
the isolated orbit using (122]) and its properties given in (1621 . fl52|) : 



1 c (0)/ \ 3 



(83) 



This is illustrated in figure [9j The "recovery distance" d from the boundary, at which 
the diverging asymptotic result goes over into the uniform result, is bounded from below 
by the Fermi wave length: 

d > h/p x ■ (84) 

For the particle number iV = 606 with p\ = 35.8186 chosen here, this corresponds to 
d > 0.03. 



+ Unfortunately, no corresponding uniform results for the kinetic energy densities have been found. 
[Note also that the relations (|23|) and (|26|) do not hold at and near the boundary r = R.] 
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Figure 10. Particle density in the two-dimensional disk billiard with N = 606, 354, 
174 and 68 particles (units: h 2 /2m = R = 1). The dotted lines are the quantum 
results, the solid lines the converged semiclassical results using the number of orbits 
and the regularizations discussed in the text. 

5. Numerical results for density oscillations 

In this section we present numerical calculations for the spatial densities in the circle 
billiard and compare our semiclassical results with exact quantum-mechanical results. In 
section I57TI we will discuss the convergence of the sum over closed orbits in semiclassical 
expressions ()14p - ffTB"]) . and in section [5721 we address the role of closed main shells and 
investigate the relation of the spatial density oscillations with the shell effects in the 
total energy of the system. 

5.1. Convergence of orbit sums 

Figure [10] shows the total particle density p(r) for four values of the number of N 
particles bound in the circular billiard. The dotted line is the quantum result ([I]), 
and the solid line the converged semiclassical result obtained as explained below. The 
numbers N = 606 and 68 correspond to filled main shells ("closed-shell systems") and 
the numbers N = 354 and 174 to half-filled main shells ("mid-shell systems"); we refer to 
the discussion of the shell structure in section [5721 (cf. figure (HJ). The agreement for the 
closed-shell systems is excellent, that for the mid-shell systems is also very satisfactory. 

In figure [TTJ we show how the convergence of the semiclassical result for N = 606 
comes about; here only the oscillating part Sp(r) of the density is shown. Dotted lines 
show the quantum results, identical in all panels, and solid lines the semiclassical results 
in various approximations. From bottom to top, an increasing number of closed orbits 
is included in the sum (|14|) . In the three lowest panels, the amplitudes Ay ([22]) for 
isolated orbits are used; in the top panel, the uniform approximations for the symmetry 
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Figure 11. Oscillating part 6p(r) of the particle density for N = 606 particles like 
in figure [TOl Dots: quantum results (identical in all panels). Solid lines: scmiclassical 
results; in the three lowest panels as isolated orbits with amplitudes (p2|) . and in the 
top panel using the uniform approximations (f66)) . (|73[) and ([78)) (see text for details). 



breaking at r = and for the pitchfork bifurcations are used. 

In the bottom panel, only the primitive radial NPOs L^ -0 ^ are included. In spite 
of their weak divergence at r = 0, their contribution already accounts for most of the 
density oscillations. In particular, the Friedel oscillations are well reproduced by the 
"+" orbit up to the highest peak; only close to the boundary its contribution diverges 
(cf. the contribution of the "+" orbit alone in figure [9]). In the second panel from below, 
the radial orbits L± -1 ^ have been added, together with the non-radial orbit A (3,1) and 
the triangle PO (3,1) emerging from by pitchfork bifurcations at r = R/3 and 

r = R/2, as shown in figure [2] (left panels). The agreement with the quantum results 
is slightly improved in the region 0.3 < r/R < 0.6, but at the cost of the divergences 
at the two bifurcation points. In the third panel, we have furthermore added the non- 
radial orbit V (2,1) and the diameter PO (2,1). This improves the agreement around 
0.1 < r/R < 0.2, but at the cost of a stronger divergence at r = 0. In the top panel, 
finally, we have used the uniform approximation (|66|) for the symmetry breaking at r = 
of the "+" and "— " orbits, including also the radial orbits L± -2 ^ and the orbits (2/c+l, k) 
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Figure 12. Same as figure [TT1 but for N = 9834 particles. 



bifurcating from them, the uniform approximation ( 1751) for the pitchfork bifurcations 
of the NPOs and POs (3,1) and (5,2), and finally the uniform approximation (173|) for 
the symmetry breaking at r = of the NPOs and POs (2k, k) including values of k up 
to k = 10. The agreement with quantum mechanics is now excellent, except very close 
to the boundary where the isolated "+" orbit (with k — 0) diverges. Employing the 
regularization at r = R discussed in section 14.31 and shown separately in figure [9] leads 
to the result shown in figure [TOl 

The same analysis is made in figure [12] for a larger system with N = 9834 particles; 
the convergence is similar as above. Here the irregular part of the oscillations is seen 
to be well reproduced by the contributions of the non-radial NPOs (2k, k) in the two 
upper panels. We have omitted the surface region r > 0.9-R because the peaks of the 
last Friedel oscillations are so high that the oscillations in the interior would not be 
recognizable on their scale. The contribution of the primitive "+" orbit alone, however, 
reproduces the quantum density almost exactly for r > 0.9R like in figure [TOl 
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Figure 13. Quantum-mechanical results (dots) and uniform semiclassical results 
(lines) for N — 9834 particles as in top panel of figure [T21 but for the oscillating 
parts of the kinetic energy densities 5£(r) (top panel), 5ti(t) (center panel), and 5r(r) 
(bottom panel). 

In figure[13]we show the oscillating parts of the kinetic energy densities of the system 
with iV = 9834 particles using the same orbits and uniform approximations as in the top 
panel of figure [T2l Like there, the agreement is best in the region < r < 0.2R, where 
the rapid oscillations are dominated by the radial orbits L± and the slow, irregular 
oscillations by the orbits (2k, k). We see that the rapid oscillations are opposite in 
phase in the two kinetic energy densities 5r(r) and Sri(r), as expressed in the relation 
( I26j) which is due to the fact that the momentum mismatch for all radial orbits has the 
constant value Q — — 1. The rapid oscillations cancel therefore in the quantity 5£(r), 
seen in the top panel of figure [131 which is due exclusively to the non-radial orbits 
(except very close to the boundary) as discussed in general in [T]. 

The small discrepancies at r > 0.2R are due to missing longer orbits in the sums 
contributing to (1141) - (fT6l) . In principle, these can be included using the general 
formulae given in section [2j There is a practical problem, however. With increasing 
values of k, the bifurcations tend to lie denser along the radial variable r. The 
uniform approximations given in section HI including those for the symmetry breaking 
at r = 0, can only be used for the regularization of a divergence around one critical 
point at a time. It is necessary that, while going away from one critical point, an 
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asymptotic domain be reached in which the uniform approximation becomes identical 
to the sum of contributions of the isolated orbits, before one approaches the next critical 
point. Midways between the two critical points, one can then switch from one uniform 
approximation to the next one. This is what we have done in the examples given 
in the figures of this section; it was possible with the limitation of /c max = 2 for the 
bifurcating orbits of type (2k + l,k). For longer orbits, successive bifurcations lie so 
close that their actions at the two bifurcation points differ by less than ~ h and hence 
the bifurcations cannot be treated separately. This is the same problem as arises for 
the pairs of orbits created by tangent bifurcations (see e.g. figure El right panels), as 
discussed in section 14.2.21 We can consider it as a lucky circumstance that, with the 
quality of the semiclassical approximation reached in the above figures, all those longer 
orbits can be practically neglected. 



5.2. Relation to shell effects in the total energy 

As already observed in figure [TUl the quality of the semiclassical approximation depends 
to some extent on the shell situation of the system with N particles. This is known also 
from semiclassical calculations of the total energy of finite fermionic systems [3T | |32 | 133] . 
To study the shell structure, it is instructive to consider the shell- correction energy 5E tot 
defined [31] as the oscillating part of the total energy: 

p\ p\ p\ 

5E tot = 2 J2 E J ~ 2 E 9( E ) dE = 2 E g(E) dE-2 Eg(E) dE . (85) 

E 3 <\ ^° ^° ^° 

Here g(E) = ^2 ■ 8{E — Ej) is the exact level density and 'g(E) its average part. A and 
A are the Fermi energies determined by N = A/"(A) = A/"(A) in terms of the "number 
counting functions" , which are defined as the first integrals of the two level densities: 

JV(E) = 2 V 9(£ - Ej) = 2 f E g(E) dE , 

3 J ° 

N{E) = 2 fg{E) dE=£-M+\ + 0{E~W) . (86) 

The explicit expression for J\f(E) is obtained from the Weyl expansion [H [9] . In figure 
IT41 we show the shell- correct ion energy as a function of the particle number N in the 
region < N < 650. A very pronounced shell structure is seen. The minima of 5E to t{N) 
correspond to closed "main shells" for which the level density at the Fermi energy has 
a minimum, too, corresponding to a low degeneracy of the spectrum. These main 
shells may be counted by the number M appearing in ( |25|) . The maxima of §E to t{N) 
correspond to regions of high degeneracy for "mid-shell" systems. This behavior of the 
energy shell structure has been studied extensively, e.g., in nuclei [35] . 

In order to calculate the shell- correction energy 5E tot semiclassically, it is convenient 
to express it directly in terms of the oscillating part 5g(E) = g(E) — g(E) of the level 
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Figure 14. Shell-correction energy SEtot(N) (|85|) versus particle number N for the 
circle billiard (energy units: Eq = 1). The circles mark the the closed-shell systems (at 
minima of <5i?tot) with N = 68 and 606, and the squares mark the mid-shell systems 
(near maxima of SE tot ) with N = 174 and 354 tested in figure ITD1 



density, for which a Gutzwiller-type trace formula [3] can be employed. SE tot (N) may, 
in fact, be approximated [311 [33] by either of the two following expressions 



SE tot (N) 



(E - A) Sg(E) dE 



5M(E) dE , 



(87) 



where 5M{E) = J\f(E) — M{E) is the oscillating part of the number counting function. 
The two relations given in (157)1 are correct up to terms of order (A — A) 2 or [5J\f(\)} 2 , 
respectively. The semi classical trace formula for 5g(E) of the circular billiard is given 
in [11]. After its integration in (I57j) . one obtains to leading order in h the following 
semiclassical approximation for the shell-correction energy in terms of the POs (v , w) 



8E tot (N) ^2El^J2i: 



fv 



sin 



PO 



7T 37T 

3v — 

2 4 



(88) 



w=1 v=2w V 2 y/7TV8in(w7r/v) 

where the lengths are given in (|47j) and f v>w = 2 — 5 V ^ W - 

In figure [15] we compare the semiclassical results (J88l) of 5E tot (N) with their exact 
values ( 1851) in the region 2075 < N < 2250. We see that the agreement is best at the 
minima (closed main shells) and fairly good around the maxima (mid-shell systems). 
In regions where 5E tot (N) has steep slopes, there are relatively large discrepancies, as 
already observed earlier [32], [33]. The reason lies in the fact that in these regions, 
the stair-case functions A(iV) and M{E) have their largest deviations from the smooth 
functions A(iV) and Af(E), so that the missing terms of order (A — A) 2 or [<5A/"(A)] 2 are 
largest, while these are smallest at closed main shells and near mid-shell regions. This 
is illustrated in figure [HH where we show the two curves M{E) and M{E) versus energy 
E in an interval covering the central part of the region of particle numbers shown in 
figure [151 They are seen to intersect at closed shells (circles) and in the mid-shell region 
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N 

Figure 15. The same as figure [TJ] in a different region of particle numbers N. Dots: 
exact values (|85p. crosses connected by the solid line: semiclassical values (|88[) . The 
circles mark closed-shell systems, the square marks mid-shell systems, and the triangles 
mark some intermediate systems tested in figure [17] below. 




Figure 16. Number counting functions Af(E) (dots connected by the solid line) and 
AT(E) (dashed line). The circles, triangles and the square correspond to those in figure 
[I5]in the region 2098 < N < 2222. 



(square). The triangles mark "difficult systems" for which the errors in the semiclassical 
results are largest, as seen in figure [T5l There is a clear correlation between these errors 
and the magnitude of 5Af(E). 
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Figure 17. Oscillating part 5p(r) of the spatial density, obtained like in the top panels 
of figures fTTI and \\% here for the six particle numbers marked by the circles, triangles 
and the square in figure \Wl 

We observe a similar trend for the quality of our semiclassical approximation to the 
spatial density oscillations. The four systems shown in figure [101 where the agreement 
with the quantum results is best, correspond to closed-shell and mid-shell systems (see 
the triangles in figure [T4l) . In figure [T71 we show the corresponding results for the five 
systems marked by circles, triangles and the square in figure [T6j Clearly, the best 
agreement is reached for the closed-shell systems (circles) . The correlation of the errors 
in 5p(r) with the magnitudes of 5M{E) is less direct than for the shell- correct ion energy, 
but the agreement in figure [UJ is clearly worst in the systems (marked by triangles) for 
which SE tot (N) has its steepest slope. The mid-shell system N = 2150 here has slightly 
larger errors (like also for 5E tot in figure [T5l) than those shown in the middle panels of 
figure [ini The same trends are also found for the kinetic energy densities. 

In conclusion we can state that, with a small number of included classical orbits 
and straight-forward uniform approximations, the semiclassical approximation for the 
spatial density oscillations in the circular billiard is excellent for closed-shell systems, 
good for mid-shell systems and still quite satisfactory for most intermediate systems. 
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6. Summary and outlook 

We have presented a complete classification of the closed periodic and non-periodic 
classical orbits in the two-dimensional circular billiard with radius R. We have divided 
the orbits into two classes: radial orbits, librating in the radial directions, and non- 
radial orbits. Analytical formulae are given for their actions, stability determinants, 
momentum mismatches and Morse indices. We have found that many of the orbits 
can undergo bifurcations at which new orbits are created. The bifurcation parameter 
hereby is their starting point r varied in the range < r < R. Periodic orbits can be 
created from non-periodic ones, and vice versa, in isochronous pitchfork bifurcations. 
Pairs of non-periodic orbits can also be generated from tangent bifurcations. Finally, a 
class of non-periodic orbits is generated from the radial periodic orbit upon breaking its 
rotational U(l) symmetry by varying its starting point from r = to r > 0. 

Employing a recently developed closed-orbit theory [Q [2], we have investigated 
the semiclassical approximation of spatial density oscillations for N non-interacting 
fermions bound in the circular billiard in terms of its closed orbits. At the critical points 
(symmetry breaking at r = 0, bifurcation points, and the boundary r = R), uniform 
approximations have to be used in order to regularize their semiclassical amplitudes. 
For the symmetry breaking and both two types of bifurcations, standard uniform 
approximations from perturbation and bifurcation theory could be implemented in a 
straightforward manner, having asserted that the bifurcations are exactly of the types 
known for periodic orbits as functions of other bifurcation parameters (such as energy, 
deformation or external field strength). Following the general analysis in [Q [2], we have 
demonstrated that also in the present system the radial orbits are responsible for the 
short-ranged regular oscillations in both particle and kinetic-energy densities, while the 
non-radial orbits create their long-ranged irregular oscillations. We have also confirmed 
here that the radial orbit with one reflection at the nearest turning point creates the 
Friedel oscillations near the boundary. 

We have finally tested the semiclassical approximation by comparing its results 
to exact quantum-mechanical densities. We investigated the convergence of the sum 
over closed orbits and the correlations of its results with the shell structure in the total 
energy. We find that, using a relatively small number of included classical orbits and 
easy analytical uniform approximations, the semiclassical approximation is excellent 
for closed-shell systems, good for mid-shell systems and still quite satisfactory for most 
intermediate systems. This is found to hold even for moderate particle numbers N > 50, 
in spite of the fact that a priori the semiclassical limit of a small Planck constant H 
corresponds to the limit of large particle numbers. 

The classification of closed orbits that we have given here is valid for spherical 
billiards in arbitrary dimensions D > 2, since the motion always takes place in a two- 
dimensional plane due to the conservation of angular momentum. The only difference 
from our present results for D > 2 is the degeneracy of the orbits: all non-radial NPOs 
starting at r > have at least one continuous symmetry of rotation and are therefore no 
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longer isolated; also the degeneracies of the POs will change (211 EE] • While the shapes 
and actions of all orbits remain the same in all dimensions, their Jacobians and Morse 
indices become different for D > 2, and the regularisation schemes discussed in section 
HJhave to be modified. This work is actually in progress. 

For smooth potentials with radial symmetry, we can expect all types of orbits found 
in the present system, except that the non-radial orbits will have smoothed reflections 
from the potential wall instead of corners. But there is reason to expect further orbits. 
There always exists a circular orbit (not present in a billiard) whose radius is given by 
the minimum of the effective potential including the centrifugal part. Further orbits 
are rather likely to be created by bifurcations, although we must expect the bifurcation 
scenarios to be different from those found here. In [37] it has, e.g., been shown that in 
radial potentials of the form V(r) = cr a , bifurcations of the circular orbit give birth to 
the (smoothed) polygon-like POs under variation of a. 

For arbitrary potentials without symmetries, it is in general rather hard to predict 
the shapes and properties of the classical orbits. In [2] we have discussed the applicability 
of our semiclassical closed-orbit theory to arbitrary potentials to some extent, and also 
given explicit results for the rectangular billiard in which all closed orbits are easily 
classified. For potentials with arbitrary shapes (smooth or hard-reflecting), one will in 
general have to rely on the numerical search of closed orbits and their bifurcations. 
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Appendix. Explicit results for specific non-radial NPOs 

All the properties of the NPOs used in the semiclassical expressions (1141) - (fT6l) for the 
density oscillations have been given in the general analytical expressions in section I3TB1 in 
terms of the angles a and (3 related by ( 13T1) and shown in figure [TJ In practice, however, 
we need them as functions of the radial variable r. Unfortunately, for the non-radial 
NPOs the inversion of the relation (I32I) is not possible in general; using the standard 
relations of the trigonometric functions this would require to find the roots of Chebychev 
polynomials. In the following, we give explicit results as functions of r for some simple 
orbits, for which the explicit inversion of ( 1321) is possible. Some of these results have 
been used in examples given in the main part of the paper. 

We recall that the properties of the radial NPOs are given explicitly in section 
13. 5} the general formulae for the Morse indices of the non-radial NPOs are given in the 
systematics of section 13.61 



Appendix A.l. v — 2, w — 1: the triangle orbit V (2,1) 

For this orbit, which exists at all starting points < r < R (see figure [8] for its shape), 
the relation (132]) becomes 

. . sin a x . , . . 

r(a) = R ; — - = R- — , x = sma. (A.l) 

v ; cos (2a) (1 - 2x 2 ) ' v ; 

Inverting the above we get 

VR 2 + 8r 2 -R / VR 2 + 8r 2 - R\ 

x(r) = , a{r) = arcsm . (A.2) 

w 4r I 4r J 

The length of the V orbit becomes 

L ^- iR ^r m i^T' (A ' 3) 



and the Jacobian becomes 



R 

J2,i(r) = — tan(2a) 
Px 



sin a + — ( 1 + 6 sin 2 a) 
R 



(A.4) 



Inserting a(r) from ( 1A.2I) into ( 1A.3j) and ( 1A.4I) . both quantities are obtained explicitly 
as functions of r. For small r we get the Taylor expansions 

L 2 ,i(r) =4fl + 2^ + 0(-^)\ J^)=^Q\oQ\ (A.5) 
Appendix A.2. v — 3, w — 1: the orbit A (3,1) 

This orbit exists only for r > r\ = R/3, where it bifurcates from the radial orbit . 
At r 

— r 3,? — Rll the triangular PO A (3,1) bifurcates from it (see the bifurcation 
scenarios in the left panels of figure [2] and the shapes of the orbit in figure [3]). At 
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the boundary r = R it becomes identical with the squared PO □ (4,1). The relation 
between r and a is given by 

r(a) = R Sma . =R- T . (A.6) 

v 7 sin (3a) (3 - 4x 2 ) v 7 

Inverting the above, we get 



1 / R 

sin Qt(r) = x(r) = - W3 , r > r\ — R/3 . (A. 7) 

The explicit expressions for its length, momentum mismatch and Jacobian are 



L 3 ,i(r) =2^p-^0R + r), (A.8) 
Qs ll (r) = ^(3r-i2)-l, (A.9) 



J3,i(r) =^(2r-i?)(3r-W— . (A.10) 
-KPa V r 

Note that, although the relations (IA.6I) and (I A. 71) only hold for r > ri, the three 

functions (1A.8H - (1A.10I) can be continued analytically to r < r x and are real in the 

entire interval < r < R. For r < r 1; we can consider them as the (real) properties 

of an imaginary "ghost orbit" ; they are required in the uniform approximation (1751) for 

the pitchfork bifurcation at r\ discussed in section 14.2. 1[ 

Appendix A. 3. v = 4, w = 1: the "pentagon" orbits P (4,1) and P' (4,1)' 

This pair of NPOs is created in a tangent bifurcation at 

r 4 ,i = 0.682489/?. (A.ll) 

Some of their shapes are shown in figure HI and their bifurcation scenario is shown in the 
right panels of figure [2j At the boundary r = R, P becomes identical to the pentagonal 
(not self-crossing) PO (5,1), and P' becomes identical to the triangular PO A (3,1). 
The connection between r and a of this orbit is given by 

r(a) = -R Sm( \ = -R-, % ^ . (A.12) 

V ; cos(4a) (l-8x 2 + 8x 4 ) V ; 

The inversion of r(x) leads to a fourth-order algebraic equation whose analytic solution 
is possible but cumbersome. Without writing it down explicitly, we can make the 
following statements. The real solutions for x(r) consist of two branches starting from 
the minimum value of r(x ) which is found at x = [(2 + ^/\Q)/12} 1 / 2 = 0.655889 and 
has the value r(xo) = r 41 given in fjA.llj) . The two branches correspond to the two 



orbits P and P' born in a tangent bifurcation at r^i. The orbit P' corresponds to 
the lower branch with x[ = 1/2 < x < x$. The orbit P corresponds to the upper 
branch with Xq < x < X\ = sin(37r/10) = 0.809017. At rf® = Rj\f2 it undergoes a 
pitchfork bifurcation, whereby the squared PO □ (4,1) is born. The lengths, momentum 
mismatches and Jacobians can in principle be evaluated explicitly but are most easily 
calculated numerically using the general formulae given in section 13.31 We note that the 



Closed orbits and spatial density oscillations in the circular billiard 



41 



analytical continuation to r < r 41 , where the orbits are imaginary "ghost orbits", here 
leads to complex quantities, which is a characteristic of the tangent bifurcation [18J. 

For the sake of an illustration, we give in the following the results for the similar 
pair of NPOs with (v,w) = (5,1), for which the two branches of x(r) can be given 
explicitly. 

Appendix A. 4- v — 5, w — 1: the pair of orbits (5,1) and (5,1)' 

These orbits are born in a tangent bifurcation at the point r 5 2 = 4/2/5. Inverting 

„ sin a —R 

r(a) = -R — ; — r = 7 = — , (A.13) 

V ; sin(5a) (5 - 20x 2 + 16x 4 ) ' V ; 

which has a minimum at xq = a/5/8, yields a bi-quadratic equation with the simple 
solutions 



sin a r = x\r) 



| y| (5 ± a/5 - 4i2/r) , (A.14) 
whereby the two signs correspond to the two orbits. Their properties are 



LM = R\J\ (3 T >/5 - 4i2/r) (4 + ^ ± 1^/5 - 4fl/r) , (A.15) 
Q 5 ,i(r) = -l + ^(5 T a/5 - 4i2/r) , (A.16) 



4i2 2 + 30r 2 - 29i2r ± r(-9i2 + lOr) a/5 - 4i2/r . (A. 17) 

The upper and lower signs in the above expressions correspond to the orbits (5,1) and 
(5,1)', respectively. At r^® = R cos(7r/5), which is a zero of the term in square brackets 
on the second line in (IA.17I) (taken with the upper sign), the NPO (5,1) gives birth to 
the regular pentagon PO (5,1) in a pitchfork bifurcation. 
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